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Abstract 

One of the main sources of error associated with the calculation of defect formation energies using 
plane-wave Density Functional Theory (DFT) is finite size error resulting from the use of relatively 
small simulation cells and periodic boundary conditions. Most widely-used methods for correcting 
this error, such as that of Makov and Payne, assume that the dielectric response of the material 
is isotropic and can be described using a scalar dielectric constant e. However, this is strictly only 
valid for cubic crystals, and cannot work in highly-anisotropic cases. Here we introduce a variation 
of the technique of extrapolation based on the Madelung potential, that allows the calculation of 
well converged dilute limit defect formation energies in non-cubic systems with highly anisotropic 
dielectric properties. As an example of the implementation of this technique we study a selection 
of defects in the ceramic oxide Li2Ti03 which is currently being considered as a lithium battery 
material and a breeder material for fusion reactors. 
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I. INTRODUCTION 



Point defects play an essential role in a number of important materials properties such 
the accommodation of nonstoichiomety and facilitation of diffusion through a crystal matrix. 
The difficulties associated with making direct observations on such small length scales mean 
it is desirable for first principles methods such as Density Functional Theory (DFT) to 
provide insight into the properties and behaviour of both intrinsic and extrinsic point defects. 

In DFT, point defects are normally modelled using the supercell methodology, whereby 
vacancy, interstitial or substitutional defects are placed in a simulation supercell which is 
then tesselated though space using periodic boundary conditions (PBCs) to create an infinite 
crystal. Therefore, any defect included in the original supercell will also be tesselated and 
the interaction of these defect images can have a significant infiuence on the defect formation 
energy. This problem is particulary acute in the case of charged defects as the Coulomb 
interaction decays slowly as a function of the separation between point charged^. A number 
of correction schemes have been devised to extract the formation energies in the desired 
dilute limit from simulations of relatively small supercells: these have been widely applied 
to systems such as silicorPEl, NaCpEI, diamoncP^, GaAPS, Inl^ and AhOp^. Inherent 
in all of these schemes is the assumption that the dielectric response of the material is 
isotropic and can be described by a single dielectric constant, e. Strictly, this only holds 
for cubic systems, but in many cases the degree of anisotropy is modest enough that the 
assumption of an isotropic dielectric response is adequatd ^ ' ^" ' ^^ '. Intuitively, one might expect 
that this would not be the case for many of the more complex crystals that are currently 
being proposed for industrial applications, particularly those with layered structures. 

One such system is lithium metatitanate, /3-Li2Ti03, which is currently under considera- 
tion for use in lithium ion batterieJ^ and as breeder material in fusion reactor^. Li2Ti03 
may be described as a distorted rocksalt structure (space group C2/c) with alternating Li, 
LiTi2 and O plane^^^^ which are clearly visible in Fig. [T] Within the LiTi2 layers the Ti 
atoms form a honeycomb structure with a Li ion at the centre of each hexagon. It is this 
layered structure that gives rise to the material's interesting dielectric properties. Currently, 
not much is known about the properties of the defects in Li2Ti03. Vijayakumar et al. de- 
termined, using empirical potentials, the relative energies required to remove the different 
Li atoms and found that the formation energy of a Li vacancy defect in the pure Li layer is 
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0.30 eV greater than in the LiTi2 layei'^. The hnear "muffin-tin" orbitals method has been 
used to study H substitution onto Li sites where the hydrogen is observed to move from the 
hthium site and bond to an oxygen forming a hydroxid^^. One of the principle reasons for 
this shortage of theoretical results is the very same problem we try to address in this paper: 
namely that the anisotropy of the system means that it is hard to extract well-converged 
formation energies. 

In this study we investigate the convergence of the formation energies of point defects 
in monoclinic /3— Li2Ti03 as a function of supercell size. Specifically, we study the V^f, 
Li^f and O^^ defects (modified Kroger- Vink notation). These defects represent a range of 
different defect types and also have high charge states and so are subject to the largest 
finite-size errors. 

II. METHODOLOGY 

The DFT simulations presented here were performed using the plane-wave pseudo- 
potential code CASTEP^°. Exchange-correlation is described using the generalised gradient 
approximation of Perdew, Burke and Ernzerhof (GGA-PBE)^^. A P-centered Monkhorst- 
PackP^ scheme was used to sample the Brillouin zone with the separation of points maintained 
as close as possible to 0.05 along each axis. 

The same pseudopotentials as in previous workpl (ultrasoft pseudo-potentials (USPs), 
generated "on-the-fly" in CASTEP, and normconserving pseudopotentials (NCPPs) from the 
standard library in Materials Studio) were employed here. The planewave kinetic energies 
were truncated at 550 eV and 1700 eV for the USP and NCPPs respectively. The Fourier 
transform grid for the electron density is larger than that of the wavefunctions by a scaling 
factor of 2.0 and the corresponding scaling for the augmentation densities was set to 2.3 
when USPs were in use. These values were determined by performing convergence tests of 
the energy from self consistent single point simulations. The lattice parameters determined 
using DFT are within 1% of the experimental as shown in Table |l} Additionally a comparison 
of all the atomic positions has been included in the Supplementary Materials. 

Supercells for the defect simulations were then constructed by creating Ixmxn repetitions 
of the relaxed unitcell in a, b and c respectively with the resulting supercells containing 
between 192 and 576 atoms. The lattice parameters and cell angles were fixed during 
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Table I. Table showing the lattice parameters and bandgaps calculated using the Ultrasoft and 
Norm-conserving pseudo-potentials compared to the available experimental data. 



Property 


Ultrasoft 


Norm-conserving 


Experimental 


Volume / 


432.98 


441.53 


427.01-12- 


a /A 


5.09 


5.12 


5.0(1^ 


h /A 


8.83 


8.90 


8.75^6 


c/A 


9.80 


9.85 


9.75^6 




100.19 


100.24 


IOO.21IS! 


Eg /eV 


3.27 


3.43 


3.9C23 



minimization of the defect containing supercells so only atom positions were relaxed until 
the residual force on each atom was < 0.08 eV A^^ and the difference in energy between 
consecutive ionic relaxation steps was < 5 x 10~^ eV atom~^. 

The calculated band gap of Li2Ti03 was 3.27 eV for the USPs and 3.43 eV for the NCPPs 
compared to an experimental value of 3.9 eV^ and values of 2.5-3.5 eV from previous first 
principles studie ^ * ^^ * ^^ l Underestimation of the bandgap is a common feature of LDA 
and GGA calculations: fortunately the lack of occupied states in the band gap for the 
fully charged defects investigated here ensures that the defect formation energies reported 
are not directly affected by this source of erroi'^, although there are still effects due to 
the localization of states. The implementation of a hybrid functional, such as the HSE 
functional^^ would most likely change the defect formation energies. However, it is important 
to note that whatever functional is used, formation energies are still subject to similar finite- 
size effects, as these are determined almost entirely by the Coulombic interaction of periodic 
images, with functional-dependent effects being limited to the polarizability and localisation 
of charges. 

Following the formalism of Zhang and NorthupP^ the formation energy of a defect is given 

by: 

Ef = ^Ifect - E^evf + ^i^^i + I^F (l) 

i 

where ^'Jefect -^Jerf the DFT total energies of a system with and without the defect, 
rii is the number of atoms added/removed, fii is the chemical potential of species i, q is the 
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charge on the defect and Ep is the Fermi energy (defined here as the valence band maximum, 
VBM). 

Representative chemical potential potentials {fiu, /^Ti and fio) were generated starting 
from the assumption that Li2Ti03 can be formed from Li20 and Ti02 via reaction [2] (this 
approach has been adopted for simplicity, though we note that this is not the traditional 
route for synthesis of Li2TiOj^, 

Li20(s) + Ti02(s) ^Li2Ti03(s). (2) 

The sum of the chemical potentials of the constituent species must equal the total Gibbs 
free energy of the Li2Ti03, i.e. 

/^Ti02 iP02,T) + /iLiaO (P02) T) = /iLiaTiOgf,) (3) 

where, /iTi02 (P02) ^) and fiu20 (P02) ^) are the chemical potentials of Li20 and Ti02 within 
lithium metatitanate as a function of oxygen partial pressure and temperature and fJ'U2TiOs(s) 
is the chemical potential of solid Li2Ti03. For a solid /i(02,T°) ^ /i(0,0) therefore 
the temperature and pressure dependencies have been dropped. Two limiting cases are 
envisaged, one in which the titanate is formed with excess Li20, ie. /XLi20 (P02) ^) = 
/iLi20(,) and /iTi02 (Poa^^) = /^Li2Ti03(,) - AiLi20(,) and similarly titania rich formation con- 
ditions where /iTi02 (P02 5 ^) = /^Ti02(3)- Here we assume Li20-rich conditions, therefore 
/iiOa (^'02'-^°) '^^^ determined from the formation energy of Li20 under standard con- 
ditions (taken from a thermochemical database"^ and the DFT total energies for Li20 and 
lithium metal. Atio2 (P02) T) can then be determined following Finnis et alW^, nu {po2)T) = 
1/2 (/^Lifo - /^|02 (P02, ^)) and similarly/iTi (P02, ^) = /^Tiol - 2/^|02 (Poi^T). For the pur- 
poses of this work a temperature of 1000 K and oxygen partial pressure of 0.2 atm were 
selected. 

In a periodic system the electrostatic energy is only finite if the total charge on the 
repeat cell is zero. Therefore, when modelling a charged defect in a simulation cell subject 
to PBCs, a uniform jellium of charge is imagined to exactly neutralize the net charge on the 
supercelpll. The electrostatic energy of a periodically repeating finite system containing a 
point charge, q, and a neutralizing background jellium is the Madelung energy, 

E = -^^ (4) 
2e ^ ^ 



where vm is the Madelung potential (for cubic systems vm = d/L, where a = 2.8373 and 
L is the supercell size length). This energy (scaled by e to represent dielectric screening in 
the material) arises due to the use of PBCs and is therefore an artifact of the simulation 
technique and must be removed from the calculated defect formation energy resulting in the 
charge correction proposed by Leslie and Gillan'^: 

ET = EAL) + (5) 

In highly ionic materials, defect charge distributions can be described as point like, so this 
correction is adequatd^, however when the defect charge distribution is more diffuse the 
correction of Makov and Payne^^ is more appropriate. Castleton and co-workers proposed 
an extrapolation procedure^^^ for the study of defects in InP whereby the limit of a fit 
to a series of formation energies obtained from supercells of increasing size was used to 
determine the formation energy of the isolated defect. However, uniformly scaling all axes 
simultaneously rapidly increases the number of atoms in the supercell. Consequently it is 
only computationally feasible to sample the smallest multiplications resulting in too few 
points to allow a reliable fit to the data. Hine et al. thus suggested an improvement to this 
scheme based on simulation of supercells comprising different multiples of the primitive cell 
along different axe^, with vm calculated separately for each cell. By plotting the defect 
formation energy as a function of vm and fitting a function of the form Ef {vm) = Ef + bvM, 
it is possible to determine the formation energy of the defect in the dilute limit from the 
intercept with the y-axis, and an effective permittivity can be extracted from the gradient 
as 6 = — g^/2e. 

The constant vm can be found using Ewald summatiorP^ 



where the sum extends over all vectors of the direct (Rj) and reciprocal (Gj) lattices, 7 
is a suitably chosen convergence parameter and is the volume of the supercell. vm is 
normally positive and hence the Madelung energy is normally negative as it is dominated 
by the interactions of the point charge and the canceling background jellium which is on 
average closer than the periodic images. For long, thin supercells this is no longer the case as 
the electrostatics are now dominated by the interactions of neighbouring point charges and 
so Vm is negative. This can, however, be viewed as an advantage since simulations can be 



performed on supercells where vm is both negative and positive and the results interpolated 
to Vm = rather than performing an extrapolation outside the range for which data is 
available. 

Fig. [2] shows the formation energy of the V^f defect as a function of Vm for a range of 
supercell shapes and sizes. The data display a wide variation and it is not possible to extract 
a single value for E^^. The origin of this variation may be deduced by examining subsets 
of the data. Shown in Fig. [2] are fits of the form Ef (vm) = + hv^ to defect formation 
energies calculated in supercells created by extrapolating in the number of repeat units along 
the h- and c-axes independently, i.e. 2xmxl (for m= 2,3 and 4) and 2xlxn (for n = 2,3 
and 4). As the effective permittivity can be related to the gradient of such a fit it is apparent 
that there is a different level of charge screening present along these crystallographic axes. 
The effective dielectric constant along b (28.3) is predicted to be more than double that 
along c (13.4). 

To account for anisotropy in the screening, the dielectric constant in Eq. [5] must be 
replaced by a tensoi'^, denoted e. For a monoclinic crystal such as Li2Ti03 the dielectric 
tensor has four non-zero components, as shown below'^ 



en ei3 

622 

ei3 £33 



(7) 



This tensor can then be incorporated into the Ewald summation to give a screened Madelung 
potential, vf;/, in the general casd^lEE 

1 erfc(7VR,.-6-i-R.) ^0 47r exp(-Gr6-G./47^) 2j 

Eq. |8] implies that it is necessary to determine the dielectric tensor for each defect cell, 
which while possible using Density Functional Perturbation Theory (DFPT)32l, is computa- 
tionally prohibitive. Here we investigate two possible methods that circumvent this problem. 

In the region far from the defect the atomic positions (and consequently the dielectric 
properties) will remain largely unaffected by the presence of the defect, however, in the region 
immediately surrounding the defect the screening properties may be strongly perturbed. If 
the perturbed region is small relative to the simulation supercell then the dielectric properties 
of the whole cell may not undergo a substantial modification. As a first approximation, we 



therefore try applying the dielectric tensor for the perfect Li2Ti03 crystal to the all defective 
systems. 

In our second approach a function Ej^vm) = — (Q'^/2)fA/ + is fitted to the defect 
formation energies determined for a number of different cell shapes and sizes. This fitting 
procedure is slightly unusual as it is the values in the x-axis that are modified by optimising 
the elements of e*^*^. Optimised values of E'^ and the associated elements of e*'^ were obtained 
using a Nelder-Mead simplex algorithm!^. 



III. RESULTS AND DISCUSSION 

The dielectric tensor for Li2Ti03 was calculated using DFPT and the norm-conserving 
pseudo-potentials. The results, presented in Table |TT) show that there is indeed a significant 
level of anisotropy in the dielectric tensor. Examining only the principal (diagonal) elements 
of e^^P""" we can see that the magnitude of e^^ is less than half that of e^*,^ and e|f . Taking the 
tensor average gives a value of 30.5, which can be compared to a value of 24 for a polycrys- 
talline Li2Ti03 sampl^^S (this value has been corrected to represent the theoretical density). 
The discrepancy between the experimental and theoretical dielectric properties may arise 
due to the inherent inability of DFT simulations to accurately reproduce experimentally 
observed band gaps. The values also deviate from those predicted by examining the subsets 
of the uncorrected data determined from Fig. [2] 

Corrections such as that of Makov-Payn^^l ^re often performed with e obtained from 
either DFPT or experiment. Fig. [s] shows defect formation energies for V^f for a selection 
of supercells as a function of f^f 5 employing e^^^P'^. The data points show that while the 
level of scatter present in Fig. |2] has been reduced there is also a poor adherence to the 
linear relationship with gradient —q^jl expected from Eq. [s] This discrepancy arises as the 
use of the dielectric tensor calculated for the perfect cell, which thus neglects the atomic 
relaxations and the consequent modification of the local screening in the vicinity of the 
defect. The modification of the dielectric properties of the supercell is further supported 
by the change in the band gap of the material upon introduction of the defect. Plotted 
in Fig. [5] are the Densities of States (DOS) for the perfect Li2Ti03 as well as the defect 
containing supercells (all DOS are produced for the 2x2x2 supercell). Fig. [s] shows that 
the bandgaps in the defect containing supercells are reduced relative to the perfect supercell 
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Table II. Effective permittivity tensor e*^ , and dilute limit defect formation energies , for several 
defect species. 



Species 


,cff 

^11 


,off 


,eff 
^33 


,eff 
^13 


Ef /eV 


LizTiOa (DFPT) 


36.1 


37.8 


17.8 


-5.0 




y-4 

Vrpj 


37.4 


37.4 


14.4 


-1.0 


6.3 




34.9 


35.0 


13.9 


-12.1 


3.5 


Of 


33.7 


55.6 


16.5 


-8.9 


7.3 



{Eg{V^t) = 2.46 eV, EgiU^^f) = 2.86 eV and ^^(O,^^) = 2.12 eV), which would suggest a 
perceptible change in the dielectric properties of the cell. 

In order to incorporate the change in the dielectric properties of the supercells induced 
by the defect, we instead fit the elements of e to give e^^. e^^ is then effectively an averaged 
picture of the dielectric tensors for the supercells included in the fit. Presented in Fig. |4] 
are plots of the formation energies as a function of fff after the fitting procedure has been 
performed for the V^f, Li^f and 0^"^ defects. The three plots in Fig. [Zjshow that by fitting 
the elements of e it is possible to substantially improve the correlation between the data and 
the relationship given in Eq. |5| thus allowing a single linear function to be fitted to the data 
and a dilute limit defect formation energy to be extracted. Residual errors associated with 
the fitting process are around 0.1 eV and likely arise either from dipole-dipole or monopole- 
quadrupole interactions not accounted in Eq. |5} or from changes in atomic configurations 
for cells with one small value of l,m,n. The relatively small errors justify our treatment of 
the defect charge state as point-like. However in complex systems where the defect charge 
is less localised, or for defect clusters, this approximation may no longer hold. 

The fitted elements of e^''^ and the resulting dilute-limit defect formation energies are 
presented in Table |llj In the case of the V^f and Li^f defects the degree of atomic relaxation 
is relatively small and the concomitant differences between (p^^^^ and e^^ are also modest. 
The level of local distortion resulting from the introduction of an O^^ defect is much greater 
than for the other defects as depicted in Fig. |6} Furthermore the reduction in the bandgap 
is greatest for this defect which is also consistent with it displaying the most significant 
perturbation in its dielectric properties. It is this higher level of relaxation that leads to 
the increased difference between e^^^PT and e^^. At larger system sizes we would expect to 
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recover values of e*' increasingly close to those of the perfect crystal, but the current results 
indicate that at the fairly small system sizes obtainable here, it is appropriate to fit to the 
observed results. 

IV. CONCLUSIONS 

In summary, we have proposed an extension of the Madelung extrapolation procedurd^^ 
for the calculation of defect formation energies in the dilute limit. This is achieved by 
incorporating the effect of charge screening, via the dielectric tensor, into the calculation of 
the Madelung potential (via Eq. [s]) and fitting the elements of the tensor and the desired 
dilute-limit formation energy to defect formation energies calculated in a range of supercells. 
We have applied the method to Li2Ti03, which has a monoclinic structure and a highly 
anisotropic dielectric tensor, and demonstrated its ability to determine defect formation 
energies converged to within around 0.1 eV even for such systems. In principle this method 
is applicable to systems of any shape and dielectric properties. Even in cubic supercells, 
local relaxation, such as that arising from defect clusters, may break the crystal symmetry 
such that the dielectric properties are anisotropic, necessitating a tensor representation of 
dielectric properties. We have also further highlighted the importance of incorporating the 
effect of lattice relaxation on the dielectric properties of the material when applying a finite- 
size correction method based on the Makov-Paynd^ approximation. 
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Figure 1. Crystal structure of /3— Li2Ti03. Yellow, green and red spheres represent titanium, 
lithium and oxygen ions respectively. The black outline represents a single unit cell. 
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Figure 2. Formation energy of the V^^^ defect for a range of supercells with differing vm- The 
wide variation in the points demonstrates that is not possible to fit a single straight line of the 
Ef (vm) = E'^ + bvM to the data. 
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Figure 3. Formation energy of Yr^f as a function of vf;f, where e^^^"'" has been used in the 
calculation of vf;f. The black dashed line represents a fit of the form given in Eq[5]to the raw data 
and the red dashed line is fitted to the corrected data. 
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Figure 5. Densities of states for perfect Li2Ti03 and the V^J^, Li^^^ and O^^ defects. The long 
dashed green lines, intermediate dashed blue lines and short dashed red lines correspond to s-, p- 
and d-derived states respectively, with the sum plotted using the solid black line. The valence band 
maximum is indicated with a dashed grey vertical line and the conduction band minimum with a 
dotted grey vertical line. 
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(a) (b) 

Figure 6. (a) Atomic arrangement surrounding an as yet unoccupied interstitial site in Li2Ti03 
and (b) same atoms after introduction of an O^'^ defect. Clearly visible in (a) is the distorted 
rocksalt structure of Li2Ti03 and the distortion arising arising due to the defect is illustrated in 
(b). It is this distortion that causes the significant change in e*'^ for this defect. All atom positions 
have been extracted from simulations in the 2x2x2 supercells. 
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